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Abstract 

We present a novel method for determining multi-fractal proper- 
ties from experimental data. It is based on maximising the likelihood 
that the given finite data set comes from a particular set of param- 
eters in a multi-parameter family of well known multi-fractals. By 
comparing characteristic correlations obtained from the original data 
with those that occur in artificially generated multi-fractals with the 
same number of data points, we expect that predicted multi-fractal 
properties are unbiased by the finiteness of the experimental data. 
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1 Introduction 

The characterisation of spatial distributions in terms of fractal concepts [0, Q 
is becoming increasingly important. In particular, many distributions in 
nature are found to have the characteristics of a multi- fractal]^, |^: among 
many examples are galaxy clustering [|l], H, strange attractors [0, fluid 
turbulence |T^, and plant distributions 
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In application, methods for estimating fractal dimensions are typically un- 
reliable through many problems. Two sources of error lie in largely unknown 
biases introduced by the finiteness of data sets, addressed by Grassberger [^], 
and in the finite range of length-scales inherent in gathered data. In many 
situations where thousands or tens of thousands of data points are known 
the biases are relatively minor; however, in some interesting problems, for 
example in the spatial clustering of underwater plants f^, only the order 
of 100 data points are known and confidence in the fractal results may be 
misplaced. 

We propose to eliminate such biases through a novel method of determin- 
ing a multi-fractal properties of a given dataset. We compare characteristics 
of inter-point distances in the data, §^ and in artificially generated multi- 
fractals. By maximising the likelihood, §|^, that the characteristics are the 
same we may model the multi-fractal nature of the data by the parameters of 
the artificial multi-fractal. By searching among, artificial multi-fractals 
with the same number of sample points as in the data, we anticipate that 
biases due to the finite sample size will be statistically the same in the data 
and in the artificial multi-fractals, §|^; hence predictions should be unbiased. 
The method also appears to give a reliable indication of the error in the 
estimates — a very desirable feature as also noted by Judd & Mees 0. 

2 Correlation density 

Borgani et al [0 have surveyed the efficacy of the most popular methods for 
determining fractal dimensions. They considered: generalised correlation in- 
tegrals, box counting algorithms, density reconstruction procedures, nearest 
neighbour methods, and minimal spanning tree methods. Their conclusion 
is that nearest neighbour and minimal spanning tree methods are inferior. 
Aesthetically this is pleasing because both of these methods discard much of 
the available information about inter-point distances that the other methods 
retain. 

Box counting algorithms are very similar to the correlation integrals ex- 
cept that they are non-isotropic. Aesthetically, an isotropic method is pre- 
ferred and so here we ignore box counting methods. 

Borgani et al also conclude that correlation integrals are most reliable 
for generalised dimensions, Dq, of positive order, g > 0, and less reliable for 
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negative order. In contrast, they maintain that the density reconstruction 
procedure is more rehable for negative order, q < 0, and less so for positive 
order. The same conclusions are also reached by Martinez Interestingly, 
both approaches are based on the same processed data; they just fit straight 
lines in complementary manners. 

Given a finite set of N data points Xi sampled from some spatial distri- 
bution, define the correlation functions Ci{r) as the fraction of points within 
a distance r of the ith data point. Then the partition function 

1 ^ 

Zir,q) = ^T.C^{rr~'^ (1) 
1=1 

is expected to scale as r"^ over the range of length-scales on which the multi- 
fractal properties hold. By fitting such a power law to Z{r, q), the generalised 
dimension Dg is then found from 

T={q-1)D,. (2) 

For example, the Hausdorff dimension is estimated with q = 0, the correlation 
dimension with q = 2, and the information dimension with g — i> 1. 

Conversely, the density reconstruction method is based on Ri{c) which 
is the radius of the smallest ball centred on the zth point that contains a 
fraction c of the data points. Then the partition function 

1 ^ 

W{c,t) = -J2R,{c)-\ (3) 
1=1 

is expected to scale as c^~'^ over a respectable range of fractions c. Then by 
fitting a power law, Dq may be again estimated from (^). 

At the heart of both of these methods are the curves Cj(r) and Ri{c) 
which are precisely the inverse function of each other. Imagine plotting all 
of the curves, i = 1, . . . ,N, in the rc-plane, and observe that: 

• the correlation integral Z{r, q) is just an average of these curves over c 
at fixed r; 

• whereas in the density reconstruction W{c, r) is an average of the curves 
over r at fixed c. 
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One method is reputably good for positive g, the other for negative q. Per- 
haps by using the data in the rc-plane directly, we can avoid such a dichotomy 
between positive and negative q. 

We propose to base our method on the correlation density, p{r, c), defined 
as being proportional to the number of Cj(r) curves (or equivalently the Ri{c) 
curves) which pass near the point (r, c) . In practise we divide the log r, log c- 
plane into rectangular bins of width Ar and Ac, and compute pjk = p{rj, Ck) 
as the fraction of curves with logc^ — Ac/2 < logCj(rj) < logc^ + Ac/2. See 
Figure |l] for example. 

To emphasise that the previous two good methods are both based on this 
same picture, observe that to O (A^): 

Z{rk,q) =Y.Pjkcr^ ^ (4) 
j 

whereas 

W{c,,T) = Y.p,krk^. (5) 

k 

3 Maximum likelihood multi-fractal 

The task now is to fit the correlation density p{r, c) in a manner so that 
subsequently we deduce a multi-fractal spectrum for the original data. Our 
idea is to determine parameters of a multiplicative multi-fractal which best 
fits the correlation density of the original data. Then our estimate of the 
multi-fractal spectrum of the data is that of the multiplicative multi-fractal. 
The best fit is determined by: 

• generating points on a random multiplicative multi-fractals with a 
given set of parameters; 

• constructing its correlation density Pjk = P{fj, Ck) on the same grid in 
the rc-plane as used in approximating p{r, c) for the original data; 

• estimating the likelihood that the two densities, p and P, come from 
the same distribution; 

• and seeking to alter the parameters of the multiplicative multi-fractal 
in order to maximise this likelihood. 
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Figure 1: perspective view of the correlation density p{r, function of 

log r and log c, for a sample of = 100 points from the artificially generated 
binary multiplicative multi-fractal shown in Figure § with parameters p = 
1/3, = 1/4. 



February 1, 1996 



Roberts Si Cronin 



6 



By probing the structure of the N data points by artificially generated fractals 
with precisely the same number of points we hope to eliminate, or at least 
reduce, biases introduced by the finiteness of the sample size. 

For definiteness, consider a sample of N data points Xj from a multi- 
fractal on the interval [0, 1] in one-dimension. 

We fit from a two parameter family of binary multiplicative multi-fractals. 
Given parameters p G [0, 0.5] and (p G [0, 0.5] a binary multiplicative multi- 
fractal is generated by the recursive procedure of dividing each interval into 
two halves, then assigning a fraction (p of the points in the interval to a 
random sub-interval of length p in the left half, and the complementary 
fraction 0' = 1 — to a random sub-interval of length p in the right half. See, 
for example, the distribution shown in Figure |^. Such a binary multiplicative 
multi-fractal has spectrum /(a) 0, §4] given parameterically in terms of 
< ^ < 1 and ^' = 1 - ^ as 

eioge + e^loge^ eiog0 + e^log0^ 

/ = i , « = i • 6 

logp logp 

Other multi-fractal properties, such as generalised dimensions, then also fol- 
low, for example, the Hausdorff dimension is just Dq = log2/log(l/p). 

The likelihood that two correlation densities come from the same distri- 
bution is determined from the statistic [0, §14.3] 

Since the likelihood is monotonic in (for fixed number of degrees of free- 
dom), maximising likelihood is equivalent to minimising x^- Of course there 
are correlations in the structure of the correlation densities p(r, c); for ex- 
ample, the Ci{r) curves are monotonic. Strictly speaking these correlations 
destroy the precise applicability of the likelihood and the statistic. How- 
ever, as is typical in estimating fractal properties, such correlations are in- 
trinsic to investigating structures over many different length scales and are 
not expected to have a serious effect. 



4 Optimum fit 



We then seek to find values of p and (p such that the comparison between 
the actual data and the artificially generated data is minimum. Unfortu- 
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Figure 2: across the middle of this graph is plotted N = 100 points on 
an artificially generated binary multiplicative multi-fractal with parameters 
p = 1/3 and = 1/4. The solid line is the cumulative mass distribution on 
the fractal, M{x) — mass in [0,x], showing regions of high density by large 
jumps in M. 
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nately, x^{P)4^) has a large amount of noise due to the random choices in 
the generation of the artificial binary multiplicative multi-fractals. Such ran- 
domness is absolutely necessary because it characterises the wide range of 
possible multi-fractals with the same multi-fractal spectrum. (It would only 
be possible to eliminate the noise if we knew a precise analytic expression for 
the correlation density P(r, c) for a binary multiplicative multi- fractal with 
parameters p and (f) and when sampled by points.) As is reasonable, the 
noise in x^(P)0) seems to be particularly prominent for small sample sizes 
N. 

The practical task is to minimise a noisy function of p and (f). Various 
approaches were tried; however, a crude procedure reminiscent of simulated 
annealing seems to be effective. We choose the m parameters {pe,(f)i), typi- 
cally m = 100, of least evaluated from a sample of parameter values. The 
sample of parameter values is initially uniformly randomly distributed over 
parameter space, [0,0.5] x [0,0.5], and then more are generated from the m 
most successful so far. Typically 20 iterations were performed, constructing 
some 2000 artificial multi-fractals and using to compare their correlation 
density with the original. Ultimately we end with a cloud of m points, as 
shown in Figure ^, corresponding to parameters which have realised multi- 
fractals with a good fit to the correlation density of the original data. The 
centre of the cloud, the mean of the parameter points, is our best estimate 
of the appropriate parameters to use in order to model the original data. 
The spread of the cloud is indicative of the sensitivity of the parameter es- 
timation to the noise inherent in fitting to the correlation density given the 
limited amount of information in data points — the spread depends upon 
the signahnoise ratio. 

The optimisation procedure reminds us of simulated annealing in that 
parameter values are retained depending upon randomness, though here the 
randomness is intrinsic to the function rather than externally imposed. The 
"temperature" , which is progressively lowered in simulated annealing, is anal- 
ogous to the threshold of the mthe best parameter value, which decreases 
from one iteration to the next. 
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Figure 3: a cloud of m = 100 parameter values which match to a good degree 
of likelihood with an artificially generated binary multiplicative multi-fractal 
with parameters p = 1/3 and 0i = 1/4. The + denotes the actual parameters 
used to generate the original data sample, whereas the o denotes the mean 
of the the cloud as our best estimate of the parameters. 
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5 Numerical experiments 

One numerical trial consists of generating a multiplicative multi-fractal with 
a specific value of its parameters, here we chose p = 1/3 and (p = 1/4 
throughout, sampled with N points where here we used N = 30, 100 (shown 
previously) or 300. This multi-fractal forms the synthetic finite data set. 
Then the above procedure is applied and the mean of the ultimate cloud of 
parameter values recorded as an estimate of the parameters of the synthetic 
data. 

To test the method, we performed 16 different trials with = 30 data 
points, = 100, and N = 300. The estimates of the multi-fractal parameters 
for = 100 and = 300 are shown in Figures ^ and ^ respectively. 
Observe that, as would be expected, the parameter predictions improve for 
the larger value of A^. Also observe that the numerical predictions of (p, (p) 
are reasonably centred upon the true values of (1/3, 1/4). The method does 
indeed seem to lack any bias due to the finite size of the sample. 

Going to an extremely small number of data points, A^ = 30, as shown in 
Figure ^ we find that apart from two bad realisations (of small p, high 0) the 
predicted multi-fractal parameters are semi-quantitative in that they indicate 
a reasonably limited area of the parameter space in which the actual multi- 
fractal parameters are to be found. We consider this an impressive result for 
such a small A^. 

A useful observation is that the spread observed between the estimates 
in such an ensemble of realisations, such as in Figure is approximately 
the spread observed in the cloud of high likelihood parameters, see a cor- 
responding cloud in Figure ^: the covariances are very similar. The similar 
spread of the cloud of high likelihood comparisons to the inherent error of the 
multi-fractal parameter prediction is also demonstrated in Figures ^ and ^ 
where the cloud for just one of the realisations is also plotted. Thus in ap- 
plication, where one generally cannot undertake more than one sample, we 
expect that the cloud will be a reliable estimate of the region within which 
we can confidentially expect the true multi-fractal parameters to lie. Indeed, 
it was almost always the case that the true value of (p, (p) lay within the 
cloud of high likelihood parameter values (the only exceptions being the two 
bad fits with A^ = 30 data points). Thus the covariance matrix of the cloud 
gives a good indication of the error in the fitted parameter values. 

We expect that the reason the cloud gives a good estimate of the error is 
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Figure 4: predicted multi-fractal parameters (p, 0), indicated by o's, from 
the maximum likelihood match to an ensemble of 16 different realisations, 
each of N — 100 data points, of a binary multiplicative multi-fractal with 
parameters p = 1/3 and = 1/4, indicated by +. The mean location of the 
realisations and predictions is indicated by a x. 
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Figure 5: predicted multi-fractal parameters {p,(l)), indicated by o's, from 
the maximum likelihood match to an emsemble of 16 different realisations, 
each oi N — 300 data points, of a binary multiplicative multi-fractal with 
parameters p = 1/3 and (f) = 1/4, indicated by +. The mean location of the 
realisations and predictions is indicated by a x. A cloud of best likelihood 
comparisons for one of the realisations is shown by the dots. Note the change 
in scale for this Figure. 
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Figure 6: predicted multi-fractal parameters (p, <;/>), indicated by o's, from 
the maximum likelihood match to an ensemble of 16 different realisations, 
each of = 30 data points, of a binary multiplicative multi-fractal with 
parameters p = 1/3 and (p = 1/4, indicated by -|-. The mean location of the 
reahsations and predictions is indicated by a x . A cloud of best hkelihood 
comparisons for one of the realisations is shown by the dots. 
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that the size of the cloud is characteristic of the signahnoise ratio in x^iPy 0)- 
However, the noise comes from randomly sampling a multi-fractal, namely 
the binary multiplicative multi-fractal, and hence is similar to the bias in- 
troduced by the same sized sample of the natural process that generated the 
original data. Thus the influence of this bias is characteristic of the influence 
of the noise we see in the analysis, and so the latter may be used to predict 
errors. 

Ultimately, researchers want to examine multi-fractal properties of the 
data. In this method these will be determined from the parameters of the 
best fit multi-fractal. For example, for any determined parameters (p, x) for 
the binary multiplicative multi-fractal, the multi-fractal spectrum of the data 
it matches would be predicted to be given by the analytic expression (|^). For 
the 16 realisations of the = 100 data-point multi-fractals, we determined 
various estimates, plotted by o's in Figure ^, of the true parameter values, 
p = 1/3 and = 1/4. For each of these estimates, we plot the corresponding 
predicted /(a) curves in Figure ^ along with the true /(a) curve. 

Observe that the predicted dimensions for low a, positive q, are quite 
good for all realisations, especially near the information dimension. However, 
predicted dimensions for high a, negative q, are poor; this is also the case for 
the Hausdorff dimension which is predicted from the maximum of the /(a) 
curve. 

6 Conclusion 

We have presented a method for predicting multi-fractal properties from ex- 
perimental data. The technique is based on maximising the likelihood that 
the given finite data set comes from a particular member of a multi-parameter 
family of well known multi-fractals. The evidence indicates that properties 
predicted by the proposed method are unbiased by the finite number of sam- 
ple points in the experimental data. 

A valuable feature of the method is that an indication of errors in the best 
fit parameters is naturally obtained. For example, the indications are that 
predictions of the generalised dimensions Dg for negative q are extremely 
unreliable. 

Perhaps the method could be improved by increasing the resolution of the 
correlation density in the (logr, logc) plane. However, this would increase the 
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Figure 7: ensemble of multi- fractal spectra /(a), dotted, for each of the 
predictions plotted in Figure ^ made from samples of = 100 data points. 
For comparison the multi-fractal spectrum for the actual fractal is plotted as 
the solid line. Observe the good estimation near the information dimension, 
but the large errors for larger a. 
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number of "bins" of low but non-zero density which may reduce the usefulness 
of the test. Furthermore, it will increase the correlations between the 
counts for each "bin" which may will reduce the quality of the statistics, 
though perhaps not significantly. 

Here we have addressed perhaps the simplest nontrivial example of the 
process. Most interesting data on spatial fractal distributions lie in two or 
three dimensions where more complex families of multiplicative multi-fractals 
need to be employed to model the characteristics of the data distribution. 
Here we only considered data in one dimension, but even then we should pro- 
pose a more adaptable approach by modelling the data through, for example, 
a family of ternary or quartcrnary multiplicative multi-fractals, for example, 
such an enlarged family of comparison multi-fractals will enable the method 
to match multi-fractals with asymmetric f{a) curves. 

Furthermore, here we knew the overall size of the fractal and that there 
is no cut-off at small scales, or indeed that there is no multi-scaling, in the 
data. The technique will have to be developed further to adapt to these 
features. 
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